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ABSTRACT 

Aperture  limited  Fourier  operators  are  like  pseudo-differential  operators,  except 
that  the  symbol  of  the  latter  is  part  of  the  unknown  kernel  function  of  the  former. 

Such  operators  arise  naturally  in  inverse  scattering  when  one  applies  an  integral  inverse 
scattering  operator  to  model  data  with  an  unknown  (such  as  a  reflection  coefficient) 
which  depends  on  the  source/receiver  location  as  well  on  interior  medium  coordinates. 
Large  wave  number  aperture  limited  Fourier  operators  can  be  analyzed  by  multi¬ 
dimensional  stationary  phase.  The  nature  of  the  asymptotic  expansion  of  such 
operators  will  then  depend  on  the  properties  of  the  kernel  function  to  which  the 
operator  is  applied.  Of  interest  in  inverse  scattering  are  distributions  with  support  at  a 
point  or  on  a  surface  (the  singular  function  of  a  surface)  and  piecewise  smooth 
functions  (whose  discontinuity  surfaces  are  reflectors). 

The  purpose  of  this  paper  is  to  present  some  features  of  the  asymptotics  of  large 
wave  number  aperture  limited  Fourier  inversion  operators  and  to  relate  these  results  to 
the  application  of  a  recently  developed  inverse  scattering  formalism  as  applied  to 
Kirchhoff-approximate  scattering  data  from  a  reflecting  surface.  It  is  shown  how  the 
latter  problem  can  be  reduced  to  the  former,  asymptotically.  Then,  the  more  simply 
derived  asymptotic  expansions  of  the  former  can  be  applied  to  predict  the  output  of  the 
latter^ 
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ABSTRACT 


Aperture  limited  Fourier  operators  are  like  pseudo-differential  operators,  except 
that  the  symbol  of  the  latter  is  part  of  the  unknown  kernel  function  of  the  former. 
Such  operators  arise  naturally  in  inverse  scattering  when  one  applies  an  integral  inverse 
scattering  operator  to  model  data  with  an  unknown  (such  as  a  reflection  coefficient) 
which  depends  on  the  source/receiver  location  as  well  on  interior  medium  coordinates. 
Large  wave  number  aperture  limited  Fourier  operators  can  be  analyzed  by  multi¬ 
dimensional  stationary  phase.  The  nature  of  the  asymptotic  expansion  of  such 
operators  will  then  depend  on  the  properties  of  the  kernel  function  to  which  the 
operator  is  applied.  Of  interest  in  inverse  scattering  are  distributions  with  support  at  a 
point  or  on  a  surface  (the  singular  function  of  a  surface)  and  piecewise  smooth 
functions  (whose  discontinuity  surfaces  are  reflectors). 

The  purpose  of  this  paper  is  to  present  some  features  of  the  asymptotics  of  large 
wave  number  aperture  limited  Fourier  inversion  operators  and  to  relate  these  results  to 
the  application  of  a  recently  developed  inverse  scattering  formalism  as  applied  to 
Kirchhoff-approximate  scattering  data  from  a  reflecting  surface.  It  is  shown  how  the 
latter  problem  can  be  reduced  to  the  former,  asymptotically.  Then,  the  more  simply 
derived  asymptotic  expansions  of  the  former  can  be  applied  to  predict  the  output  of  the 
latter. 
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This  is  another  in  a  series  of  papers  [Bleistein,  1987a,b]  which  analyzes  by  classical 
methods  an  important  advance  in  inverse  scattering  proposed  by  Beylkin  [1985].  The 
method  applies  to  a  variety  of  source/receiver  configurations  in  a  fairly  arbitrary  host 
medium  with  the  viability  of  the  experimental  arrangement  characterized  by  the 
nonvanishing  of  a  certain  Jacobi  determinant.  Beylkin  confirms  the  validity  of  his 
technique  through  the  application  of  pseudo-differential  operator  theory,  which  predicts 
the  validity  of  the  result  in  a  large  wave  number  limit.  Thus,  those  features  of  the 
unknown  function  which  dominate  the  large  wave  number  Fourier  transform  of  the 
function,  will  dominate  the  output  when  his  inversion  is  applied  to  model  data.  For 
example,  when  the  unknown  is  piecewise  smooth,  the  discontinuities  of  the  unknown 
function  (reflectors)  will  dominate,  followed  by  discontinuities  of  the  first  derivative, 
etc. 

Bleistein  [1987a],  proposes  a  modification  of  the  Beylkin  inversion  operator 
designed  to  process  data  from  piecewise  smooth  functions  and  produce  a  reflector  map. 
That  is,  at  each  discontinuity  of  the  medium  parameters,  this  operator  produces  a 
band-limited  delta  function  which  peaks  on  the  discontinuity  surface  —  the  singular 
function  of  the  surface  —  scaled  by  a  factor  from  which  the  change  in  medium 
parameters  across  the  surface  can  be  estimated.  That  paper  deals  with  the  case  in 
which  only  the  propagation  speed  varies  and  Bleistein  [1987b]  addresses  the  problem  in 
which  both  propagation  speed  and  density  vary  across  the  reflector.  Current  research 
focuses  on  isotropic  and  anisotropic  elastic  wave  field  inversion.  The  inversion  operator 
takes  the  form  of  an  integration  over  the  data  collection  surface.  It  has  the  structure  of 
a  Kirchhoff  migration  operator  [Schneider,  1978]  and,  therefore,  is  called  Kirchhoff 
inversion. 

In  both  of  the  above  cited  papers  by  this  author,  the  validity  of  the  inversion 
operator  was  confirmed  by  applying  the  operator  to  Kirchhoff-approximate  data  from  a 
single  reflector  and  analyzing  the  output  by  multi-dimensional  stationary  phase  applied 
simultaneously  to  the  integration  over  the  data  collection  surface  and  the  reflecting 
surface.  Stationarity  conditions  from  the  latter  integral  merely  impose  Snell’s  law  on 
the  ray  trajectories  from  a  source  and  receiver  to  the  reflecting  surface.  Stationarity 
conditions  on  the  former  integral  ties  the  integration  variables  to  the  output  point  in  a 
manner  that  is  dependent  on  the  source/receiver  configuration. 

The  main  point  of  this  paper  is  that  the  analysis  of  this  multi-fold  integral  can  be 
reduced  by  classical  methods  to  the  analysis  of  large  wave  number  aperture  limited 
Fourier  transform-like  integrals  of  functions  which  depend  on  both  the  spatial  and  wave 
vector  variables. 

This  is  exactly  the  structure  of  a  pseudo-differential  operator.  However,  in 
application  to  inverse  scattering  (in  particular,  when  the  modeling  data  is  Kirchhoff 
data),  the  symbol  of  the  pseudo-differential  operator  is  part  of  the  unknown  function. 
Thus,  important  results  from  that  theory  do  not  seem  to  apply  here.  Furthermore,  the 
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classical  analysis  of  these  integrals  leads  to  detailed  information  about  the  structure  of 
the  output  of  this  type  of  inversion  for  different  classes  of  functions  of  interest  in 
inverse  scattering.  In  particular,  the  relationship  between  the  peak  amplitude  of  the 
output  and  reflectivity  is  accessible  in  the  classical  theory  and  has  not  yet  been 
deomonst rated  in  the  pseudo-differential  operator  theory. 

Aperture  limited  Fourier  transforms  have  the  advantage  that  the  phase  of  the 
integrand  is  linear  in  both  the  spatial  and  Fourier  variables,  making  the  application  of 
the  method  of  stationary  phase  much  easier  than  it  is  for  the  analysis  of  the  inverse 
scattering  operator.  In  the  latter,  the  phase  is  a  sum  of  differences  of  travel  time 
trajectories.  The  sum  is  from  source  and  receiver  to  either  reflection  surface  point  or 
output  point.  The  trajectories  are  the  rays  of  geometrical  optics  characterized  in  the 
general  case  only  by  the  differential  equations  which  they  satisfy. 

The  following  section  presents  the  asymptotic  analysis  of  Fourier  integral 
operators  on  piecewise  smooth  functions  of  the  spatial  variable,  but  with  the  additional 
feature  that  they  depend  on  the  unit  vector  in  the  dual  transform  variable,  as  well. 

This  is  followed  by  some  numerical  examples  which  demonstrate  the  theory. 

The  fourth  section  introduces  the  application  of  Kirchhoff  inversion  to  Kirchhoff 
data  from  a  single  reflector.  It  is  shown  there  that  consideration  of  output  points 
anywhere  but  in  a  neighborhood  of  the  reflector  can  be  neglected  under  reasonable 
assumptions  and  then  that  the  integration  over  the  reflector  need  only  be  carried  out 
“nearby”  this  output  point.  In  this  small  region,  there  exists  a  transformation  of 
coordinates  from  frequency  and  spatial  variables  over  the  observation  surface  to  wave 
vector  variables.  This  transformation  has  the  property  that  it  transforms  the  multi¬ 
traveltime  phase  function  of  the  integral  operator  into  the  Fourier  dot  product  of  wave 
vector  with  spatial  variables.  At  that  point,  the  integral  is  reduced  to  one  of  the  type 
considered  in  the  previous  sections  and  the  asymptotic  output  is  known  from  the 
results  of  those  sections. 

Thus,  this  type  of  travel  time  inversion  is  seen  to  be  a  generalized  forward  and 
inverse  Fourier  transform  of  a  function  of  both  spatial  and  wave  vector  variables.  This 
is  suggested  in  Beylkin  [1985],  although  that  author  seems  to  prefer  a  space-time 
formulation  of  the  solution  and  interpretation  of  the  result  as  a  generalized  Radon 
transform. 

APERTURE  LIMITED  FOURIER  IDENTITY  OPERATORS 

Consider  the  integral 

1  =  /  dZk  I  dZx'  f(x'/L,k)  exp  [»fc*(x-  *)1  .  (1) 

(2tt)3  jD,  L  J 

A 

Here,  x,  x'  and  k  are  three  component  vectors  and  £  is  a  unit  vector.  If  the  domain  Z?* 


were  infinite  in  extent,  and  the  function  /  were  independent  of  k,  then  the  integral,  /, 
would  just  be  equal  to  /  ( x/L )  for  x  in  Dz-  and  zero,  outside,  for  a  broad  class  of 
functions  /.  This  is  just  Fourier  inversion.  The  structure  of  the  analysis,^  then,  comes 
from  the  nature  of  the  domain,  Dk,  and  the  additional  dependence  of  /  on  k. 

The  function  /  is  assumed  to  have  as  many  derivatives  as  necessary  to  carry  out 
all  the  differentiations  below.  For  functions  that  are  only  piecewise  smooth  enough, 
decompose  the  domain  of  integration  into  separate  domains  whose  boundaries  include 
all  of  the  discontinuities  of  /  in  x\  Then,  the  integral  over  each  subdomain  is  of  the 
type  defined  here.  It  will  be  seen  below  that,  asymptotically,  the  integral  depends 
primarily  on  the  boundary  values  of  the  integrand  at  certain  critical  points.  The 
manner  in  which  the  critical  points  are  determined  will  make  it  clear  that  for  a  sum  of 
such  integrals,  the  output  will  depend  on  the  jump  in  the  integrand  across  these 
discontinuity  surfaces.  That  is,  a  critical  point  on  a  boundary  surface  for  one  integral 
will  simultaneously  be  a  critical  point  for  the  other  integral  sharing  the  same  piece  of 
boundary. 

The  length  scale,  L,  is  assumed  to  characterize  the  size  of  the  derivatives  of  /. 
That  is,  derivatives  with  respect  to  y'  =  x'/L  should  be  comparable  in  size  to  /,  itself. 
For  convenience,  the  same  parameter  L  is  used  to  characterize  the  length  scales  of  the 
bounding  surface  of  Dx-.  For  example,  L  might  be  a  “typical”  principal  radius  of 
curvature  or  a  lower  bound  of  the  principal  radii  of  curvature  for  the  boundary  of  Dx. 

In  application  to  inverse  scattering,  the  domain,  Dk,  is  symmetric  with  respect  to 
the  origin.  That  is,  whenever  k  is  in  £>*,  so  is  -  k.  The  reason  for  this  is  that  k  is 
proportional  to  a  frequency,  cj,  thus  including  —  k  in  Dk  whenever  it  includes  +  Jfc. 
Thus,  symmetry  of  D *  will  be  assumed  below,  although  it  is  not  essential  to  the 
analysis.  Furthermore,  it  is  assumed  that  Dk  does  not  contain  the  origin.  In  fact, 
denote  by  K  the  minimum  distance  from  the  origin  to  the  domain  Dk.  Now  introduce 
the  dimensionless  variables, 


p  =  k/K ,  y'  =  x'jL ,  X  =  KL  , 


(2) 


and  rewrite  (1)  as 


/(A,»)  = 


_A_ 

27T 


f  p2  dp  sin0  do  d<j>  f  d3y'  f(y',  p)  exp  i'A  pp-(y-y') 
JDf  JD,-  L 


(3) 


In  this  equation,  p,  6,  <j>,  are  polar  coordinates,  and  the  unit  vector,  p,  is  given  in  terms 
of  0  and  <t>  by 


p  =  (sinflcos^,  sin0sin<£,  cos#)  . 


(4) 
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The  domains,  Dy-  and  Dp,  are  the  images  of  Dx-  and  D*,  respectively,  under  the  scaling 

(2) .  In  particular,  the  minimum  distance  from  the  origin  in  p  to  Dp  is  equal  to  unity. 

The  objective  now  is  to  analyze  /{ A,y)  asymptotically  for  large  A  as  y  varies  over 
some  domain  which  includes  Dy-.  This  Is  aperture- limited  Fourier  inversion  when  the 
aperture  is  such  that  the  length  scales  of  /  and  its  support  in  the  x-domain  are  “many” 
wavelengths  for  all  of  the  available  information  in  the  k-domain.  Formally,  this  is  the 
asymptotic  expansions  of  /  ( A,y)  as  A  — +  oo. 

The  product,  Ap  >  A  >  >1,  appears  as  a  “natural”  large  parameter  in  the  integral 

(3) .  Thus,  multi-dimensional  stationary  phase  in  the  five  variables,  y',  6,  <f>  will  be 
carried  out,  with  Ap  as  a  formal  large  parameter,  followed  by  integration  with  respect 
to  p. 

The  Significance  of  the  Boundary  Values  in  Dy- 
Define 

$  =  p.(y_y')  (5) 

and  note  that 

V'$  =  -p*=0  .  (6) 


Here,  V'  is  the  gradient  of  $  with  respect  to  the  three  variables,  y\  Since  this  gradient 
is  never  zero,  there  can  be  no  stationary  points  of  the  fivefold  integral.  Thus,  proceed 
to  calculate  the  asymptotic  expansion  of  the  integral  (3)  by  using  integration  by  parts 
(the  divergence  theorem).  To  do  so,  first  rewrite  the  integrand  as 

I 

»Ap 


/(y\  p)  exp 


*Ap$ 


'V'' 


[p  /(y'.p) 


exp 


iAp$ 


+  /i(y'»P)  exp 


»Ap$ 


/i(y\p) 


(7) 


When  this  is  substituted  into  (3)  and  the  divergence  theorem  is  applied  to  the  first 
term,  the  result  is 


/(A,y)  =  /0(A,  y)  +  — -  /i(A,  y)  , 


(8) 
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'»(*.  *>  -  T3F 


I  p  dp  sinO  dd  d(j>  •  I  dSy- h-p  f(y‘,p)  ex p  «Ap$  , 

J  D„  Js,  L  J 


/j  (A,  y)  =  —  f  pdpsmddOd<t>  f  d3y'  f\{y,p)  exp  t‘Ap$ 

2 *  J  d,  J  Dt 


In  this  equation,  n  denotes  the  normal  to  the  boundary,  Sy-,  of  the  domain  Dy-.  It 
should  be  noted  that  the  integral  /1(A,y)  is  exactly  like  7(A,y),  itself.  However,  in  (8), 
it  is  multiplied  by  1/i'A.  Now  consider  the  effect  of  repeating  this  process  with  /i(A,y). 
We  would  again  obtain  an  integral  like  /o(A,y),  but  multiplied  by  another  power  of 
1/iA.  Continuing  recursively,  we  obtain  an  asymptotic  series  of  integrals  over  Dp  and 
Sy-.  Thus,  it  is  reasonable  to  conclude  that  the  leading  order  asymptotic  expansion  of  7 
must  come  from  the  analysis  of  /o(A,y),  unless,  of  course,  /  (y',p)  were  identically  zero 
on  Sy- .  In  fact,  this  last  observation  leads  to  the  following  result. 

Lemma  1: 

Suppose  that  /  (y',p)  is  infinitely  differentiable  in  Dy-  and  vanishes  “infinitely 
smoothly”  on  the  boundary,  Sy-.  Then  7(A,y)  is  asymptotically  zero  to  all  orders 
of  1/A. 


Proof: 


In  (8),  70(A,y)  is  zero  because  of  the  assumptions  on  /.  However,  the  assumptions 
on  /  are  true  for  /lt  as  well.  Thus,  repeat  the  integration  by  parts  process  and 
obtain  another  integral  like  7(A,y)  but  multiplied  now  by  1/(«A)2,  with  an 
integrand  f?,  which  also  satisfies  the  conditions  placed  on  /.  Repeat  the  process 
recursively  and  obtain  any  desired  algebraic  power  of  1/iA  as  a  multiplier.  This 
completes  the  proof. 


An  Example 


As  an  example  of  this  result,  consider 


0,  z'3  <  0, 


/(*•)  = 


fL  /i'3  exp  -  L  fx3  ,  *3  >  0  . 


Note  that  this  function  is  infinitely  differentiable  for  all  x .  Assume  that  Dz-  is  all 
space  and  that  D *  is  the  symmetric  domain  K  <  \  k3  \  <  K±.  After  integrating  in  all  of 
the  variables  except  fc3,  7,  as  defined  by  (1),  is  given  by 
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-  K 


1  =  ~^=  l  r  +  /, 


L 

Vn  J  -  AT, 


K, 


dk , 


exp 


-  ^e,’r/4^n**  +  ,'^32  +  iV/4s5nA:3 


(11) 


In  this  equation, 


(12) 


This  result  can  be  obtained  by  using  results  about  the  modified  Bessel  function  to  be 
found  in  Watson  [1980]  and  Abramowitz  and  Stegun,  [1965].  From  this  result,  it 
follows  that 


/=o 


exp 


-  Vn 


A  =  KL 


(13) 


This  estimate  is  independent  of  z.  With  a  little  more  effort,  one  can  show  that,  for 
example,  when  A  =  3,  |  / 1  <  .04,  independent  of  z.  On  the  other  hand,  directly  from 
(10),  /  (0,0,7, /5)«1.83.  Thus,  the  aperture-limited  Fourier  inversion  provides  a  poor 
approximation  of  this  infinitely  differentiable  function  for  “large”  values  of  A. 

Remark: 

This  example  demonstrates  the  lemma  (although  the  infinite  domain  of  the 
example  requires  a  little  extra  effort).  In  this  example,  the  decay  is  actually 
exponential  in  V\.  In  general,  one  can  only  predict  “faster  than  algebraic”  decay. 

The  point  of  the  lemma  is  that  the  integrand  in  (8)  only  has  boundary  critical 
points  and  then,  only  if  the  function  /  (y\p)  does  not  vanish  infinitely  smoothly.  That 
is,  it  is  the  discontinuities  of  /  that  dominate  the  integration.  (A  boundary  point  where 
/  does  not  vanish  infinitely  smoothly  is,  after  all,  a  discontinuity  of  the  function  or  one 
of  its  derivatives,  since  its  interior  limit  is  not  equal  to  its  exterior  limit.)  Thus,  for  this 
class  of  functions,  the  large  wave  number  aperture-limited  Fourier  integral  operator  in 
(l)  is  not  an  identity  operator,  but,  at  best,  it  is  an  operator  whose  output  depends 
only  on  the  discontinuities  of  /  (or  its  derivatives)  in  some  as  yet  undetermined  way. 

Stationary  Phase  Analysis  for  I0 

Consider,  now,  the  asymptotic  analysis  of  70(A,v),  in  (9).  Again  the  phase  is 
given  by  (5),  except  that  y  is  a  function  of  two  suri.xe  parameters,  say,  a !  and  cr2. 
Thus,  fourfold  stationary  phase  in  the  variables,  olt  o2,  0,  and  <f>  will  be  applied  to  this 
integral.  The  first  derivatives  of  $  are  given  by 
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Ooi 


»  =  1,2, 


d<f> 

dO 


*  #9<f> 

If 


sin0£-(y-y')  , 


9  =  (cos#cosd>,  cos#sin<£,  —  sintf)  , 


£  =  (—  sin<£,  cos<£,  0)  . 


(14) 


Note  that  the  vectors  0  and  $  are  orthogonal  to  p.  Setting  these  four  first  derivatives 
equal  to  zero  has  the  following  geometrical  interpretation.  At  the  stationary  point,  p 
must  be  orthogonal  to  two  linearly  independent  tangent  vectors  in  Sy-.  Thus,  p  and  n 
must  be  colinear  or  anti-colinear.  Furthermore,  y  —  y  is  orthogonal  to  two  linearly 
independent  tangent  vectors  on  the  unit  sphere  in  p  at  the  stationary  point.  Therefore, 
p  and  y  —  y  are  colinear  or  anti-colinear.  Consequently,  p,  n,  and  y  —  y  must  all  line 
up. 

Given  y,  drop  a  perpendicular  to  5y-.  This  determines  a  point,  y',  on  Sy-  and  a 
corresponding  pair,  al  and  a2.  Furthermore,  for  this  point,  h  and  y  —  y  line  up.  Now 
find  p  to  line  up  with  these  two  vectors.  This  determines  a  choice  of  0  and  <t>.  In  this 
manner,  Oj,  o2,  9,  and  <f>  are  all  determined  as  functions  of  y.  Given  the  stationary 
four-tuple  (al,a2,&,4>)^  a  second  set,  (alfa2,n  —  6,ir  —  0),  also  satisfies  the  stationarity 
conditions,  merely  replacing  the  vector  p  in  the  stationary  set  of  vectors  by  —  p.  When 
p  is  a  direction  in  Dp,  —  p  is,  as  well.  Thus,  there  may  be  none,  one,  or  more  than  one 
such  pairs  of  stationary  points  for  a  given  choice  of  y.  Those  values  of  y  for  which 
there  are  no  stationary  points  are  points  at  which  J0(A,y)  (and,  therefore,  7(A,y),  as 
well)  is  asymptotically  of  lower  order  than  for  those  points  where  there  are  stationary 
pairs  of  four-tuples.  Stationary  points  exist  when  the  aperture  in  p  contains  the  normal 
from  y  to  one  or  more  points  on  Sy-. 

The  method  of  stationary  phase  requires  the  computation  of  the  determinant  and 
signature  of  a  certain  matrix,  the  Hessian  of  the  phase  function.  For  simplicity  of 
notation,  introduce 


0  —  03  >  $  —  » 


(15) 


and  then  define  the  matrix 


- 

d2$ 

<9er,  do. 

*\  J 


1  -  4  . 


(16) 


Here,  all  derivatives  are  to  be  evaluated  at  the  stationary  point.  This  calculation  is 
carried  out  in  Appendix  A.  The  result  is 
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»)  ~  S  °(»>  *’  w J (*• »)  • 


(17) 


In  this  equation,  the  summation  is  over  stationary  points  and  y',  r*  and  p  are  to  be 
evaluated  at  each  stationary  point.  The  function  J(A,y)  is  given  by 


i(A,v)  =  ^ 


f  ie. 

J  t'p 


exp 


*ApM3  I  if  -  y'l  +  */W4 


(18) 


with  p  and  p3  defined  in  equation  (20),  below.  The  domain  of  integration  is  the 
portions  of  the  rays  in  the  stationary  directions,  ±p,  which  lie  in  Dp.  The  function  G  is 
given  by 


G(y,  y',p) 


f(y,p) 

|  l  -  PiPzKi  |y-y')  |  |  1  -  p2p3K2  |y-y'|  | 

77 F 

(19) 


Here,  j  —  1,2,  are  the  principal  curvature  vectors  on  Sy-  at  the  stationary  point  and 
K],  j  =  1,2,  are  their  magnitudes.  Furthermore, 


Pi  =  sgn  p  ■  ,  p2  =  sgn  P  ‘  k2  ,  =  sgnp  ■  (y  -  y')  ,  p  =  sgn 


(20) 


In  carrying  out  the  sum  over  stationary  points,  note  that  when  p  is  replaced  by 
-  p,  p3  changes  sign,  so  that  the  exponent  changes  sign.  Furthermore,  the  factor 
n  p/p  changes  sign,  so  that  the  quotient,  n-p/p  changes  sign.  Thus,  these  two 
integrals  add  up  to  yield  an  integral  with  a  signed  variable  p  on  both  positive  and 
negative  parts  of  the  real  line.  Below,  let  us  assume  that  this  pairing  has  occurred. 
The  integration  is  to  be  carried  out  over  the  two  intervals  simultaneously  and  the 
summation  is  over  only  half  of  the  stationary  points,  say  the  ones  for  which  p3  >  0. 
When  p3  =0,  some  equivalent  ordering  in  pj  or  p2  must  be  used. 


y  on  Sy- 


Consider  the  integral  in  (18)  when  y  approaches  the  surface  Sy-.  Then,  for  the 
nearest  stationary  point,  the  phase  approaches  zero  in  this  limit.  (There  may  be  other 
stationary  points  further  away  on  Sy-  for  which  the  limit  is  not  zero.)  For  this 
distinguished  stationary  point,  J(A,y)  is  0(1)  in  A  in  this  limit,  whereas  it  is  0(1/A) 
for  all  other  stationary  points.  Thus,  J(A,y)  changes  its  order  in  A  when  y  moves  to 
Sy-.  Therefore,  it  is  reasonable  to  expect  the  sum  to  be  dominated  by  this  nearest 
stationary  point. 
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It  is  shown  in  Appendix  A  that  for  y  near  Sy-,  p  at  this  distinguished  stationary 
point  is  equal  to  zero.  Then  J  can  be  recognized  as  a  band-limited  step  function.  The 
magnitude  of  the  step  can  be  determined  by  evaluating  G(y,y',p)  in  the  limit  when  y 
is  on  Sv ■ .  The  result  is 


Gpe.k  =/(»,?)>»  on  S„- 


(21) 


That  is,  the  magnitude  of  the  step  is  just  the  value  of  /  (y,p),  which  is  the  magnitude 
of  the  discontinuity  of  /  across  the  surface  Sy'  at  the  stationary  value  of  p.  For  two 
integrals  over  domains  sharing  the  same  segment  of  boundary  surface,  the  stationary 
point  is  shared  as  well.  In  fact,  the  only  differences  in  evaluation  of  the  integral  arise 
from  the  different  values  of  the  amplitude  function,  G  (y,  y',p)  in  the  two  integrands 
and  the  fact  that  n  has  opposite  direction  in  the  two  integrals.  Thus,  the  sum  of  the 
integrals  will  yield  a  difference  of  function  values  which  reduces  to  just  the  jump  in  the 
function  /  when  y  is  actually  on  the  surface  Sy’  and  for  the  same  value  of  p.  This  is  the 
main  result  regarding  the  asymptotic  identity  operator,  (3). 

It  should  be  noted  that  if  the  domain  Dp  does  not  contain  the  stationary  value 
p  =  n,  that  is,  if  the  normal  to  Sy-  at  a  particular  y  is  not  a  direction  in  Dp,  then  there 
is  no  stationary  point.  In  this  case,  the  asymptotic  expansion  of  /  is  lower  order  and, 
presumably,  smaller  in  magnitude.  This  will  be  demonstrated  in  the  next  section 
through  numerical  examples. 

Example:  Aperture  Limited  Fourier  Inversion  of  a  Step  Function 


As  a  simple  example  of  this  type,  consider  the  function 


A,  1  <  y'z  <  2, 
0,  otherwise. 


(22) 


Assume  that  Dy ■  is  the  support  of  /  and  that  Dp  is  a  domain  symmetric  with  respect  to 
the  origin  that  intersects  the  line  Pi  —  Pi  =0  with  p  >  1  on  Dp.  Call  that  intersection 
D3.  By  substituting  this  function  into  the  definition  (3)  of  /( A, y) ,  it  is  fairly 
straightforward  to  obtain  the  result 


-  £  l 


D,  *P  3 


exp 


*'Ap3  (jfi  - 1) 


-  exp 


(»3  -2) 


(23) 


The  integral  here  can  be  recognized  as  a  difference  of  band-limited  step  functions  which 
define  the  support  of  the  function  /  (y',p)  in  the  y ’-domain. 
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For  this  example,  /i,  defined  in  (20),  is  identically  zero,  as  are  the  principal 
curvatures  at  both  boundaries.  Furthermore,  n-p  =  ±l  at  the  two  boundaries, 
accounting  for  the  difference  in  sign  in  the  two  terms  of  the  integrand.  Finally,  note 
that  the  effect  of  integration  in  y‘,  followed  by  px  and  p2  is  to  evaluate  p  =  (0,0, ±1), 
that  is,  colinear  or  anti-colinear  with  the  normals  to  boundary  of  the  support  of 
/  (y',p).  Thus,  if  we  were  to  introduce  a  multiplier,  g(p)  into  the  definition  of  /,  all  of 
the  integrations  could  still  be  carried  out  to  yield  (23),  except  that  now  the  integrand 
would  have  an  additional  factor  of  g  (0,0,1)  as  predicted  by  the  theory.  In  summary, 
the  asymptotic  result  is  exact  for  this  simple  example. 

Extracting  Information  about  /  on  Sy- 

In  applications  of  this  analysis  to  inverse  problems,  the  aperture  limited 
information  about  /  (x'/L,  k)  in  the  Fourier  domain  is  the  known  data  and  the 
objective  is  to  extract  information  about  the  function  /  ( x/L ,  k)  in  the  spatial  domain. 
The  analysis  presented  here  suggests  that  from  large  wave  number  aperture-limited 
data,  one  can,  at  best,  expect  to  extract  information  about  the  boundary  values  of 
/  [x/L,k)  at  a  value  of  k  distinguished  by  the  fact  that  it  is  colinear  with  the  normal  to 
the  boundary  passing  through  x.  More  generally,  one  might  hope  to  determine  the 
value  of  the  jump  in  /  ( x/L,k )  across  its  discontinuity  surfaces  in  the  x-domain  at  the 
distinguished  k  value. 

There  are  two  issues  to  be  addressed:  first,  how  one  might  most  easily  extract  that 
information  about  /.  In  practice,  band  limited  step  functions  are  not  easy  to  recognize. 
That  is,  given  numerical  output  rather  than  an  analytic  expression  for  the  band-limited 
step,  the  location  of  the  step  and  its  magnitude  are  not  easily  extracted,  nor  is 
information  about  the  normal  to  the  discontinuity  surface  of  /.  In  fact,  the  amplitude 
of  the  band-limited  step  function  is  actually  equal  to  zero  right  at  its  discontinuity, 
significantly  mitigating  the  effect  of  growth  of  the  integral  /  from  0(  1/A)  to  0(1)  in 
this  region.  Thus,  searching  for  the  midpoint  of  a  band-limited  step  would  not  seem  to 
be  the  most  desirable  approach.  The  second  issue  to  be  addressed  is  how  one  could  go 
about  determining  the  distinguished  value  of  p  =  k  without  having  to  actually  construct 
the  discontinuity  surface.  These  issues  will  be  addressed  in  the  order  presented. 
However,  note  that  neither  of  the  solutions  to  these  problems  is  new.  The  former 
problem  was  first  discussed  in  Bleistein  and  Cohen,  [1979]  and  the  latter  was  addressed 
in  Bleistein,  [1987a, bj. 

Processing  for  a  Scaled  Singular  Function  of  the  Boundary  Surface,  Sy- 

The  singular  function  of  a  surface,  to  be  denoted  by  7(x),  is  a  Dirac  delta  function 
with  support  on  the  given  surface.  Thus,  it  is  a  distribution  having  the  property 
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(24) 


Here,  S7  is  the  surface  having  7(2)  as  its  singular  function.  Knowledge  of  7(1) 
constitutes  mathematical  imaging  of  the  surface  S1.  A  graphical  display  of  7(x) 
provides  an  actual  image  of  the  surface. 

With  minor  modification  of  the  kernel,  the  integral  operator  in  (1)  can  be 
transformed  into  one  which  the  singular  function  of  the  boundary  surface  as  its  output, 
within  a  scale  factor  which  will  provide  information  about  the  boundary  values  of  the 
function  /(y,p). 

What  we  need  to  transform  the  integral  J(A,y)  in  (18)  into  a  band-limited  delta 
function  is  a  multiplier  which  cancels  the  factor  h-p/ip  of  the  integrand.  Thus,  we 
need  a  multiplier  of  ip  on  half  of  the  p-domain  and  —  ip  on  the  image  of  that  domain 

A  . 

through  the  origin.  The  multiplier,  iksgnii'k  with  u  a  constant1  unit  vector,  will  do 
the  trick.  When  this  factor  is  inserted  into  the  integrand  in  (1),  it  will  have  exactly  the 
desired  effect.  After  rescaling, 

ik  sgn  u  •  k  =  t  Ap  sgn  u  •  p  .  (25) 


A 

When  Dk  contains  all  directions,  k,  this  multiplier  is  not  defined  in  the  plane  of 
directions  through  the  origin  with  tt  as  normal.  In  most  applications,  however,  there  is 
at  least  one  plane  of  directions  in  which  there  is  no  information.  That  is,  the  domain 
Dk  will  exclude  some  plane  through  the  origin.  Choose  u  as  the  unit  normal  to  that 
plane. 

Proceeding  with  this  idea,  in  place  of  the  integral  (l),  consider  the  following: 

I  = — i—r  f  ksgnukd3k  f  d3xf(x'/L)ex p  ik-(x-x’)  (26) 

(2tt)3  JDk  JD,  L  J 


‘There  is  no  reason  why  U  must  be  a  constant  vector.  I  simply  have  not  found  a  need  to  make  tt 
variable. 


-  II  - 


After  rescaling  according  to  (2)  and  (4) ,  one  obtains  in  place  of  (3)  the  integral, 


I  (A,  y)  =  - ~  f  p2  sgn  u  •  p  dp  sin#  dB  d<f> 

(2*) 

•/  d3y'f{y',p)  exp  [»App*(y-y') 

*  n 


and  the  integral  I0  in  (8)  and  (9)  is  replaced  by 


/  o(A,»)  = 


2ir 


I  p2  sgn  u  ■  p  dp  suiO  dO  d<f> 

JD, 


(27) 


dSy  h  p  /(f\  p)  exp 


tAp$ 


(28) 


The  asymptotic  analysis  of  this  integral  proceeds  as  for  /o,  yielding  in  place  of  (17)  and 

(18), 


r(x,V)~5]G(,,lr%«7(A,,) , 


-  r  r 

J  =  ~2tT  J  dpexp  ,/isAp  +  *W4 


pA  =  n  -p  sgn  u  ■  p 


(29) 

(30) 


Except  for  the  new  parameter,  (±1),  the  constituent  fimctions  and  parameters  here 
are  still  defined  by  (19)  and  (20). 

For  y  on  Sv-,  there  is  a_stationary  point  at  y'  =  y  as  long  as  n  is  a  direction  in  D^. 
For  this  stationary  point,  J  =  0(A),  whereas  for  all  other  stationary  points,  or  for  y 
not  on  S„',  J  =0(1)  in  A.  Thus,  as  y  approaches  Sy-,  one  stationary  point  dominates 
the  value  of  /( A,y)  and  the  function  value  is  larger  by  0(A)  than  the  value  for  y 
bounded  away  from  5y-.  As  noted  earlier,  for  this  distinguished  stationary  point,  p  =  0, 
at  least  for  y_near  Sv- .  Recalling  that  we  must  sum  over  the  two  stationary  points,  ±p, 
the  integral  J  becomes 


J  (A,y) 


Aa*4  r  p~ 

2*  J-P+ 


i\p  |  y  —  y ' 


(31) 
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The  limits  of  integration,  p_  and  p+  are  the  intersections  of  the  ray  from  the  origin  in 
the  stationary  direction  p  with  the  domain,  Dp.  In  the  k-domain,  that  is,  undoing  the 
scaling  defined  by  (23), 


J  (A,  y) 


2tt 


dk  exp 


(32) 


Here,  k±  have  definitions  completely  analogous  to  p±.  In  either  form,  J  can  be 
recognized  as  a  symmetric  (zero  phase)  band  limited  delta  function  with  support  at 
|  y  —  y#  |  =0  [equivalently,  |  *  —  x  |  =0],  which  occurs  when  y  is  on  Sy-  [or  x  is  on  an 
equivalent  surface  in  the  x-domain,  5X  ].  In  this  limit,  the  dominant  term  in  (29)  is 
readily  evaluated  with  the  aid  of  (21): 


/(A,  y)  ~  A(p+—  p_)/x4/(y,p)/?r  =  {k  +  -  k_)nAf{x/L,  k)/irt  y  on  Sy  .  (33) 


Thus,  the  value  of  /  on  Sy-  is  proportional  to  the  interval  width  in  the  k-domain  along 
an  appropriate  ray,  multiplied  by  p4/(x/£,i)/j r.  If  the  integrand  contained  a  filter 
factor,  F[k)  =  F(Kp),  then  the  factor  p+  —  p_  is  replaced  by  the  area  under  the  Filter 
function  in  the  stationary  direction  of  p. 

Qualitatively,  for  x  near  Sx',  the  dominant  term  in  the  sum  in  (29)  has  the  form  of 
an  aperture-limited  singular  function  scaled  by  a  slowly  varying  function.  The  scale 
factor  becomes  the  jump  in  the  function  /  multiplied  by  when  x  is  on  Sx-  and  the 
aperture-limited  singular  function  becomes  the  area  under  the  band  pass  filter  in  a 
distinguished  direction,  divided  by  2ir. 

Within  a  sign,  this  output  is  the  leading  order  high  frequency  band  limited 
asymptotic  expansion  of  the  normal  derivative  of  the  discontinuous  function,  /  (y,  p)  in 
the  x-domain.  We  obtain  this  result  even  though  we  use  no  a  priori  knowledge  about  the 
boundary  surface  or  its  normal  direction  in  the  processing  operator,  I . 


Example:  Processing  for  Singular  Function  Output 


We  return  to  the  example  (22).  However,  now,  instead  of  applying  the  operator 
(l),  we  apply  (26).  At  the  same  level  of  computation  as  in  (23),  the  output  is  now 


-  „  A  A  sgn  u3  r 

■ 

1 

1  (A,  y)  -  J  dp 3 

2  7T  J  D, 

exp 

»Ap3(y3-i) 

-  exp 

*Ap3(y3-2) 
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A  sgnuz 


sin 


-Wy-i)  sin 


^P3(y3  -2) 


y3-i 


y3-2 


P3+ 


(54) 


P  3  —  P  3- 


One  can  see  here  that  /  (A,  y)  is  a  difference  of  band  limited  delta  functions  that  peak 
on  the  boundary  surfaces  of  Dy-.  That  is,  the  output  is  proportional  to  the  band 
limited  singular  functions  of  the  boundary  surface(s),  Sy-.  The  scaling  factor  is  the 
jump  in  /,  namely,  A  multiplied  by  ±nz/ic.  The  choice  of  sign  here  comes  from  the 
fact  that  the  outward  normal  to  the  boundary  has  opposite  direction  on  the  two 
portions  of  the  boundary  surface. 

For  A p3_  “large”  say,  at  least  3,  the  peaks  of  the  two  singlular  functions  will  be 
well  separated.  If,  in  addition,  the  percentage  bandwidth,  (p3+  —  P3-)/(p3+  +  P3-),  is 
“large,”  say  bigger  than  .5,  the  side  lobes  of  the  sine  functions  appearing  in  (34)  will  be 
well  damped  and  the  boundaries  will  be  well  delineated  by  the  band  limited  singular 
functions.  The  peak  values  of  /  are  given  by 


/  (A,i1,x2,1) 


A  A  sgn  u3  /  v  Asgnuz  /§  t  s 

(P3+  —  P3-)  —  (*3+  “  *3-) 


I  (A,Xj,z3,2) 


A  A  sgn  u3 
n 


(p3+  ~Pz~) 


A  sgn  u3 

7T 


(*3+  ~*3-) 


(35) 


Thus,  the  modified  operator,  (26),  produces  a  result  from  which  the  boundary  surfaces 
and  the  amplitude  of  the  discontiuity  of  /  in  the  x-domain  is  more  readily  determined 
than  from  the  ordinary  Fourier  inversion  in  (1). 

It  should  be  noted  again  that  if  the  aperture  Dp  [equivalently,  Z?*]  did  not  contain 
some  segment  of  the  line  pl  =  p2  =0,  then  /(A,  y)  =0.  Since  this  line  is  the  direction 
of  the  normal  to  the  discontinuity  surface  of  /,  this  result  is  consistent  with  the  theory. 


Example:  Aperture  Limited  Fourier  Inversion  for  a  Ramp 


The  next  example  will  show  two  features  of  aperture-limited  large  wave  number 
Fourier  inversion  that  are  predicted  by  our  results.  Consider  the  function 


/(y'.p) 


'A  (1  —  y'3),  0  <  y'3  <  1 
0  ,  otherwise  . 


(36) 


with  Dy-  the  support  of  /( f',p). 
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It  is  straightforward  to  carry  out  all  of  the  integrations  of  the  operator  I  defined 
by  (1)  except  p3  to  obtain  the  result 


'<A-  »>-■£/ 


dpz 


2*  jd,  *P  3 


exp 


«^P3  y3 


+ 


1 

»'Ap3 


exp 


*Ap3(y3  - 1) 


-  exp 


*Ap3  y3 


(37) 


If  D3  were  the  entire  line,  then  the  Fourier  transform  of  the  first  exponential  would  be 
a  step  with  support,  y3  >  0.  The  next  two  terms,  which  correct  this  step  to  produce 
the  appropriate  finite  ramp,  are  lower  order,  0(  1/A),  for  large  wave  number  aperture- 
limited  data.  Note  that  a  phase  depending  on  the  difference,  y3  —  1,  arises  only  in  this 
lower  order  term.  The  reason  for  this  is  that  the  function  defined  by  (36)  is  continuous 
at  this  boundary. 

From  aperture-limited  large  wave  number  data,  one  should  not  expect  to  detect 
the  lower  order,  0(1/A),  contribution  to  the  Fourier  transform  of  /  (y,  p),  but  only  the 
leading  order  term.  To  leading  order,  then,  7(A,y)  is  a  band  limited  step  function  with 
discontinuity  on  the  discontinuity  surface  of  /  (y,  p)  and  amplitude  equal  to  the 
amplitude  of  the  discontinuity  of  /  (y,  p).  We  could  use  the  operator,  7  (A,  y),  defined 
by  (26)  or  (27),  to  more  easily  detect  the  location  of  this  discontinuity  surface  and  the 
value,  A.  In  this  case,  the  leading  order  step  would  be  replaced  by  a  leading  order 
singular  function  of  the  boundary  surface,  y3  =  0,  while  the  second  order  ramp  would 
be  replaced  by  a  second  order  step  on  the  interval,  0  <  y3  <  1. 

Of  course,  one  might  consider  using  yet  another  operator,  with  multiplier  -A  2p2 
to  detect  the  jump  in  the  first  derivative  of  /  (y,p).  While  this  will  work  for  synthetic 
examples,  it  would  be  considerably  less  reliable  with  field  data  and  “many” 
discontinuity  surfaces  of  /  (y,p)  and  its  derivatives. 


Determination  of  p  at  the  Distinguished  Stationary  Point 


We  will  show  here  how  to  determine  the  distinguished  value  of  p  at  the  stationary 
point  by  direct  computation  of  an  operator  output.  We  consider  this  method  superior 
to  determination  of  the  distinguished  p  graphically  from  the  image  of  the  boundary 
surface.  This  method  has  been  used  successfully  by  Parsons  [1986]  and  Sullivan  and 
Cohen  [1987], 

The  determination  of  p  is  equivalent  to  the  determination  of  sin#  or  cos#  and  sin^ 
or  cos <f>  at  the  stationary  point.  To  fmd  these  angles,  note  first  from  (4)  that 
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COS0  = 


^3 

k 


El 

p 


cos4>  =  —  1 

y/  k\  +  k$ 


Pi 

\/Pl  +  P2 


(38) 


Suppose  that  we  define  two_new  integral  operators  with  these  factors  introduced  into 
the  kernels  of  the  operator,  I .  That  is,  starting  from  (26),  define 


Is  —  f  k3sgnu-kd3k  f  d3x  f{x/L,  k)  exp  ik-(x-x’) 
JDt  JD,  L 


=/, 


kkisgnk 3  ,  f  , 

-  d3k  f  d3x  f  (*' /L,  Jk)exp 


ik  •  (1  —  1') 


(39) 


(40) 


Let  us  now  consider  the  asymptotic  analysis  of  and  /^.  Clearly,  it  will 
proceed  exactly  as  above,  with  only  the  amplitudes  of  the  results  being  affected  by  the 
changes  introduced  here.  Thus,  the  peak  values  of  these  new  integral  operators  can 
also  be  predicted.  They  will  differ  from  the  result  for  /,  itself,  by  the  factors, 
cos 0,  cos <f>,  each  evaluated  at  the  distinguished  stationary  point.  That  is, 

cos„  =  !iz*  ,  co8*  =  lises. ,  y  on  S/  .  (41) 

l  peak  l  peak 

The  function,  sin^,  can  be  determined  this  way  as  well,  while  sin#  can  be  defined  in 
terms  of  cos 0,  taking  the  positive  square  root  (0  <  0  <  7r).  Given  these  values  at  the 
stationary  point,  p  is  determined  there. 

Integrands  with  Other  Types  of  Singularities 


It  should  be  noted  that  the  multiplier,  ik  sgn  u-k,  is  a  “best  choice”  only  because 
the  properties  of  the  function,  /  ( x'/L ,  k),  namely,  that  this  function  was  assumed  to 
be  smooth  but  not  to  vanish  smoothly  on  the  boundary  of  its  domain  of  definition.  In 
particular,  the  original  operator,  /,  will  optimally  depict  a  function  which  is  a  (sum  of) 
Dirac  delta  function(s).  Furthermore,  I(x)  will  also  be  a  better  operator  than  I  for 
singular  functions.  That  can  readily  be  verified  now,  as  follows.  Let  us  consider  the 
application  of  /  as  defined  by  (1)  to  the  function, 


fix’ /L,k)  =  7(z' /L)  , 


(42) 


with  7 (y'),  the  singular  function  of  a  surface,  S,  as  defined  in  equation  (24)  and  the 
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related  discussion.  Substitution  of  this  function  into  (1)  and  use  of  the  change  of 
variables  defined  by  (2)  and  (4)  leads  to  the  result 


/(A,y)  = 


2?r 


f  p2dpsin0  dd  d<j>  f  dSex p 
J  D,  J  S 


t'Ap$ 


(43) 


with  $  defined  by  (5). 

Thus,  without  introducing  an  additional  multiplier  to  modify  the  operator  /  to  / , 
this  integral  is  exactly  like  the  integral  I  0(A,y)  defined  by  (28),  with  the  amplitude  of 
that  integral,  sgn(u-p)  h’p  f  (y',p)=fj.if  (y\p)  replaced  by  unity.  Consequently,  the 
asymptotic  analysis  for  this  integral  has  already  been  done.  The  output  will  be  an 
aperture  limited  version  of  ^(y).  That  is,  it  will  be  a  band  limited  delta  function  along 
each  normal  which  is  a  direction  from  the  origin  to  D *.  The  structure  of  the  delta 
function  will  depend  on  the  extent  of  that  ray  in  the  k-domain  and  also  on  any 
smoothing  that  might  be  applied  through  the  introduction  of  a  high  wave  number  band 
pass  filter  along  that  ray. 

Asymptotically,  then,  the  high  wave  number  aperture  limited  Fourier  inversion  of 
data  for  a  singular  function  behaves  just  as  the  exact  result  did  for  a  planar  surface. 

Minor  Extension 


In  the  applications  of  interest,  the  function,  /,  may  depend  on  x,  in  addition  to  the 
other  dependencies  indicated  in  (1).  This  does  not  change  the  results  stated  here,  since 
x  merely  acts  as  a  parameter  with  respect  to  the  integrations  in  (1). 

A  Form  That  Arises  in  the  Analysis  of  Kirehhoff  Data  for  the  Inverse  Scattering  Formalism 


A  last  extension  to  be  considered  here  is  the  case  in  which  the  integrand  has  both 
a  singular  function  and  a^  smoother  amplitude  function.  Namely,  in  (1),  replace 
/  ( x'/L ,  k)  by  /  (x'/L,x/L,  k)~j(z'/L).  Then,  after  exploiting  the  singular  function,  the 
integral  (1)  becomes 


/=rV  I  d*k  I  dS'f(x/L,x/L,/k)exp\ik-(x-x) 
(2?r)3  JDt  Js, 


(44) 


or,  in  dimensionless  variables, 
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/  = 


dSy-  f(y‘,  y,  p)  exp 


»Ap$ 


(45) 


Here,  $  is  defined  by  (5) .  Except  for  the  fact  that  we  have  not  used  polar  coordinates, 
the  integral  is  of  the  same  form  as  (28)  with  the  amplitude 
n-p  sgnu-pf  (y‘ ,p)  =  (y\p)  replaced  by  /(y',y,p).  The  analysis  proceeds  as 

before,  with  the  extra  dependence  on  y  in  /  merely  playing  the  role  of  a  parameter  as 
regards  the  asymptotic  analysis  of  the  integral. 

Consequently,  the  asymptotic  expansion  of  the  integral  in  (45)  is  given  by  (29), 
with  the  following  insertions.  The  function  G (y, y',p)  is  defined  by  (19),  except  that 
/  (y',p)  is  replaced  by  /(y',y,p)  and  J  is  given  by  (30),  (31)  or  (36)  with  replaced 
by  unity.  In  particular,  for  x  near  5,  the  asymptotic  expansion  is  dominated  by  a 
scaled  band  limited  singular  function.  When  z  is  on  S,  the  scale  factor  is  just  the 
amplitude,  f  (z’/L,x/L,k),  evaluated  at  the  stationary  point  and  determinable  as 
described  in  the  discussion  above  for  determination  of  p. 


Summary 


It  has  been  shown  here  that  the  high  wave  number  aperture  limited  Fourier 
inversion  of  the  Fourier  transform  of  a  piecewise  smooth  function  is  dominated  by  the 
function  values  on  the  discontinuity  surface(s).  The  inversion  is  approximately  a  band 
limited  step  function  of  normal  distance  from  each  point  on  the  discontinuity  surface. 
The  amplitude  of  the  step  function  is  proportional  to  the  jump  in  the  function  across 
the  surface  at  the  point  in  question.  By  modifying  the  inversion  operator,  the  output 
can  be  transformed  into  the  singular  function(s)  of  the  discontinuity  surface(s)  of  the 
original  function,  again  scaled  by  the  jump  in  the  original  function  at  each  point  of  the 
discontinuity  surface.  This  modification  facilitates  the  numerical  identification  of  the 
discontinuity  surface(s)  and  the  amplitude  of  the  jump  at  each  point  on  the  surface(s). 

This  aperture  limited  inversion  of  data  was  shown  to  produce  an  aperture  limited 
approximation  of  the  singular  function. 

In  the  language  of  asymptotic  expansions  of  integrals,  given  an  output  point,  y, 
the  high  wave  number  aperture  limited  Fourier  inversion  of  a  piecewise  smooth 
function  or  of  the  singular  function  of  a  surface  is  dominated  by  certain  critical  points 
which  can  be  determined  as  functions  of  y.  A  critical  point  consists  of  a  location  y  in 
the  spatial  domain  and  a  direction  p  in  the  dual  Fourier  domain.  When  y  is  on  a 
discontinuity  surface  of  the  smooth  function  or  on  the  support  surface  of  the  singular 
function  and  the  normal  to  the  surface  in  question  is  a  direction  p  in  the  aperture  Dp, 
the  output  is  at  least  0(A)  larger  than  it  is  otherwise.  This  increase  in  order  can  be 
exploited  to  detect  the  discontinuity  surface(s)  of  piecewise  smooth  functions  or 
support  surface(s)  of  singular  functions  and  estimate  an  amplitude  function  on  the 
support  surface(s). 
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NUMERICAL  EXAMPLES 


We  show  here  some  numerical  examples  that  demonstrate  the  theory.  All  are  two 
dimensional,  allowing  us  to  use  the  third  dimension  to  depict  the  amplitude  of  the 
output. 

Aperture  Limited  Singular  Function 


As  a  first  example,  consider  the  singular  function  of  the  line,  xl  =  0.  This  function 
has  as  exact  Fourier  transform,  <5(Ar2).  Thus,  the  aperture  limited  Fourier  inversion  will 
be  nonzero  only  if  D *  contains  some  segment  of  the  line  k2  =  0.  In  the  context  of  our 
theory,  the  normal  to  the  support  domain  of  the  singular  function  is  (0,1).  Therefore, 
this  must  be  a  direction  in  Dk.  When  the  converse  is  true,  the  theory  predicts  that  the 
output  will  be  asymptotically  zero  to  leading  order.  In  this  simple  example,  the 
asymptotic  result  is  exactly  true. 

Figure  1  depicts  the  singular  function  of  the  line,  =  0.  This  depiction  of  the 
singular  function  extends  over  only  one  sample  point  and  is  therefore  a  full  bandwidth 
singular  function  for  the  discrete  (FFT)  Fourier  transform  to  be  used  below.  The  extra 
“spike”  at  the  comer  is  shown  for  normalization  purposes;  it  will  appear  in  subsequent 
figures.  Figure  2  shows  the  k-domain,  Dk,  to  be  used  for  inversion.  Note  that  this 
domain  includes  the  normal  direction  to  the  line,  X\  =0  of  Figure  1.  Figure  3  shows 
the  aperture  limited  inversion.  The  aperture  limited  singular  function  adequately 
identifies  the  original  function.  Figure  4  depicts  another  k-domain,  Dk,  for  inversion  of 
the  same  data.  Note  that  this  domain  does  not  include  the  normal  direction  to  the  line, 
z1  =0,  although,  on  a  percentage  basis,  it  is  much  larger  than  the  previous  domain. 
Figure  5  is  a  plot  of  the  aperture-limited  inversion  for  this  domain,  Dk,  on  the  same 
scale  as  Figures  1  and  3.  Figure  6  shows  the  same  output,  scaled  up  by  six  orders  of 
magnitude. 

The  more  general  case  is  demonstrated  in  Figures  7-11.  Figure  7  depicts  the 
singular  function  of  a  circle.  Figure  8  depicts  a  k-domain,  Dk,  in  which  the  angle  of  k  is 
not  restricted,  but  only  its  magnitude.  Figure  9  is  the  inversion,  showing  an  aperture 
limited  Dirac  delta  function  along  every  radial  line  (every  normal  to  the  original  curve). 
Figure  10  shows  an  alternative  k-domain,  Dk,  for  which  the  direction  of  k  [equivalently, 
k j  is  restricted,  as  well.  Figure  11  shows  the  inversion.  Here,  the  singular  function  is 
adequately  defined  only  for  the  normal  to  the  circle  in  the  angular  aperture  of  Dk  and 
the  output  falls  off  strongly  to  zero  outside  of  that  angular  aperture. 
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Relevance  to  Inverse  Scattering 


This  last  example  has  relevance  to  inverse  problems.  The  surface,  S7,  might  be  a 
reflector.  Thus,  determination  of  x/L )  from  aperture-limited  Fourier  data 
constitutes  mathematical  imaging  of  the  reflecting  surface.  In  practice,  the  singular 
function  of  the  surface  might  also  be  multiplied  by  a  reflection  coefficient  which  is  a 
function  of  k.  As  noted  above,  k  coligns  with  the  normal  to  57  and  therefore  it  is  fixed, 
as  well.  Because  of  the  subtleties  of  the  dependence  of  the  reflection  coefficient  on  k  in 
applications,  the  reflection  coefficient  that  is  determined  need  not  be  the  normal 
reflection  coefficient.  See  Bleistein  (1987a,  1987b]  and  the  discussion  of  the  next 
section. 


APERTURE  LIMITED  INVERSE  SCATTERING 


From  references  cited  here,  as  well  as  many  others,  the  evidence  is  overwhelming 
that  high-frequency  inverse  scattering  is  closely  related  to  aperture-limited  Fourier 
inversion.  However,  it  is  not  always  true  that  high-frequency  inversion  is  equivalent  to 
large  wave  number  inversion.  See,  for  example,  Beylkin,  Oristaglio  and  Miller  [1985]; 
Miller,  Oristaglio  and  Beylkin  [1987];  and  Mora  [1987].  For  simple  scattering 
experiments  in  which  the  sources  and  receivers  are  on  the  “same  side”  of  the  scattering 
region,  high  frequency  does,  indeed,  correspond  to  large  wave  number.  That  is  the 
type  of  experiment  of  interest  here. 

The  inversion  operator  in  equation  (49),  below,  was  proposed  by  the  author 
[Bleistein,  1987b]  for  an  inverse  scattering  problem  with  the  following  properties: 

(1)  a  variable  background  sound  speed  and  density  for  back  propagation  are 
prescribed; 

(2)  a  source  receiver  array,  such  as  (a)  common  (or  single)  source,  multi-receiver, 
(b)  common  receiver,  multi-source,  (c)  common  offset  between  source  and 
receiver,  is  prescribed; 

(3)  the  data  can  be  described  as  high  frequency  data. 

Under  further  conditions  to  be  discussed  below,  the  operator  was  shown  to  yield  a 
reflector  map  and  a  means  of  estimating  parameter  changes  across  the  reflector.  In  this 
section,  we  show  how  analysis  of  the  application  of  the  inversion  operator  to  model  data 
is  related  to  the  discussion  of  large  wave  number  aperture  limited  Fourier  inversion  as 
discussed  in  the  previous  two  sections. 
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The  Inverse  Problem 


Consider  an  inverse  scattering  problem  for  acoustic  wave  propagation  in  which  the 
source/receiver  pairs,  x,  and  xr,  respectively,  are  identified  by  a  parameter  ^=(^1,^2) 
as  follows: 


xs  —  *«(£)>  xr  —  xr(£)  • 


(46) 


For  example,  in  the  case  of  a  common  (or  single)  source  experiment,  xs  would  be  a 
constant  vector  denoting  that  fixed  position  and  the  function  xr(£)  would  be  a 
parametric  representation  of  the  receiver  surface.  For  the  common  receiver  case,  the 
roles  of  xs  and  xr  would  be  reversed.  For  the  experiment  with  common  (or  fixed)  offset 
(between  source  and  receiver),  both  xs  and  xr  would  vary  with  £. 

It  is  assumed  that  the  sound  speed  c  (x)  and  the  density  p(x)  are  known  down  to  a 
reflecting  surface  5;  below  5,  these  variables  take  on  new  values,  e  +  (x)  and  p+(x),  to 
be  determined,  along  with  the  surface  S,  itself. 

It  is  further  assumed  that  tt  (x,x4,w)  is  the  response  to  an  impulsive  point  source 
at  xs.  Above  5,  u  satisfies  the  wave  equation 


pV  • 


(47) 


Below  S,  u  satisfies  this  equation  with  right  side  zero  and  with  p  and  c  replaced  by  p+ 
and  c  +  ,  respectively.  Furthermore,  u  must  satisfy  two  continuity  conditions  on  S, 
namely  that  u  and  1  /pdu/dn  are  continuous  across  S. 

The  total  solution  above  S  may  be  viewed  as  the  sum  of  a  free-space  Green’s 
function  and  the  upward  scattered  response  to  the  reflector.  The  Green’s  function 
satisfies  (47)  everywhere,  with  c(x)  analytically  continued  throughout  space  as  an 
infinitely  smooth  function.2  Denote  the  second  term  by  us(x,xt,u).  Thus,  the  data  for 
the  inverse  scattering  experiment  is 


D[£,u)  =  u5 


(48) 


In  Bleistein  [1987b],  the  following  inversion  operator  to  process  this  data  was 


2Clearly,  there  are  extensions  in  which  the  reference  is  only  “piecewise  smooth”  and  the  Green’s 
function  is  replaced  by  a  “primaries  only”  Green’s  function,  or  one  that  contains  multiple  reflections 
through  interfaces  above  S. 
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iu>  dui  F(w)  exp 


r(z,  i,)  +  r(z,  ir) 


D((,u)  . 


(49) 


The  domain  of  integration  is  the  set  of  £- values  which  are  required  to  cover  the 
source/receiver  array.  The  domain  of  integration  in  w  is  limited  by  the  high  pass  filter 
F(w).  The  functions  r(x,x,)  and  A(x,xt )  (r(x,xr)  and  A  (x,xr)]  are  the  WKBJ  or  ray- 
theoretic  phase  and  amplitude  of  the  constant  density  Green’s  function  with  source  at 
x ,  [ xr ]  and  observation  point  x. 

The  function,  h(x,£),  is  the  essential  element  of  this  inversion  formula.  It  is  the 
determinant, 


V  r(x,  X,)  +  r(x,  xT) 

a  V 
dZl 

t(x,  x9)  +  r(x,  a 

r)] 

-*-V 

dt2 

r(x,  x,)  +  r(x,  I,) 

(50) 


The  factor,  Nl(x,^),  was  omitted  in  the  above  cited  references.  This  is  a 
neutralizer  or  cutoff  function  having  the  property  that  h  is  of  one  sign  on  the  support  of 
this  function.  The  neutralizer  is  equal  to  unity  except  near  the  boundary  of  its  support 
domain,  where  it  decreases  C°°  smoothly  to  zero.  The  introduction  of  this  function 
precludes  spurious  artifacts  in  the  output  arising  from  the  boundary  of  S^.  Further 
discussion  of  the  assumption,  h  =  0,  is  provided  in  Beylkin  [1985]  and  in  Beylkin, 
Oristaglio  and  Miller  [1985] . 

The  derivation  of  this  inversion  formula  starts  from  Beylkin’s  [1985]  equation 
(4.3),  expressing  that  result  in  the  frequency  domain.  Asymptotic  analysis  indicates 
that  for  high  frequency  band  limited  data,  Beylkin’s  inversion  formula  will  produce  a 
band-limited  step  function  wherever  the  impedance  has  a  discontinuity.  That  is,  the 
output  of  the  inversion  formula  is  similar  to  the  output  of  aperture  limited  Fourier 
inversion  applied  to  piecewise  smooth  functions.  Proceeding  with  an  _analysis 
equivalent  to  the  determination  of  the  modified  Fourier  inversion  operator,  /,  of  the 
previous  sections,  leads  to  a  modification  of  the  Beylkin  inversion  operator  that  yields  a 
singular  function  output  instead  of  step  function  output.  The  key  to  this  modification 
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is  the  identification  of  an  appropriate  wave  vector,  k,  which  is  also  indicated  in 
Beylkin’s  [1985]  paper. 

In  Bleistein  [1987a, b],  the  validity  of  the  modified  inversion  was  confirmed  by 
applying  the  operator  to  model  data  derived  using  the  Kirchhoff  approximation.  The 
purpose  of  this  section  is  to  show  that,  in  fact,  the  analysis  can  be  reduced  to  the 
simpler  problem  of  aperture  limited  Fourier  inversion  as  discussed  in  the  previous 
sections. 

To  show  how  this  is  done,  the  inversion  formula  is  applied  to  the  upward  scattered 
response  from  a  single  reflector,  S,  modeled  by  the  Kirchhoff  approximation.  This 
representation  can  be  found  in  many  sources,  including  Bleistein  [1986],  eq.  (74).  In  the 
notation  used  here  the  result  is 


tOJ 


P(xr) 

pM 


1/2 


J  R{x',xs)A(x',xa)A(z\xr)h-V'  r(x',  x,)  +  r(x',  xr) 


•  N2(x‘)exp 


r(z4,  x)  +  r(xr,  x) 


dS' 


(51) 


In  this  equation,  V'  denotes  a  gradient  with  respect  to  the  x  variables  and  R  (x',x,)  is 
the  geometrical  optics  reflection  coefficient, 


R{x',*s)  = 


1  Mx\*s) 

p(x')  dn 


1 

1  .  1  ..  , 

dr(x',x,) 

2  ‘ 

p+ix‘) 

4(x')  C2(z') 

dn 

1 

dr(z',zt) 

+  1 

1  .  1  .  . 

dr(x‘,x$) 

2  ‘ 

P{x') 

dn 

P+(x‘) 

400  c2(x') 

dn 

.  (52) 


The  unit  normal  n  points  upward  and  d/dn  =  n-V';  N2{x)  is  a  neutralizer  function, 
equal  to  unity  except  in  a  neighborhood  of  the  boundary  of  S.  If  S  is  of  infinite  extent, 
let  N2(z’)  vanish  “sufficiently  far”  from  the  origin.  The  reason  for  this,  is  that,  again, 
edge  effects  from  the  boundary  of  the  domain  of  integration  are  not  of  interest  in  the 
current  analysis.  The  result  (51)  is  substituted  into  (49)  to  obtain  the  following  multi¬ 
fold  integral  representation  of  the  output  (3(x )  when  applied  to  this  synthetic  data: 
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0{x)  ~  -j?  I di( 


h(x,t)\N  i(x,e) 


A(z,  zt)  A(z,zr)  Vr(x,  xg)  +  Vr(x,  X,) 


J  w2  du;  F(u/) 


/  R{z\  Zg'j  A(x  \  zt )  A(x',  Zr)exp  j*w#(*  z\z„  zr ) 


JV. 


,(*')  n  -V'  |r(x',xg)  +  r(z',zr) 


dS’  . 


(53) 


In  this  equation, 


$(x,x',xs,xr)  =  r(x',xg)  +  r(x',xr)  -  r(x',xg)  +  r(x',zr) 


(54) 


is  the  difference  of  travel  times,  source  point  [xg]  to  input  point  [x']  to  receiver  point 
[xr]  minus  source  point  to  output  point  (x]to  receiver  point.  The  surface  S  is  described 
parametrically  in  terms  of  two  parameters,  Oj ,  ct2,  by 

x  =  z‘(a),  o  =  [p |,*a)  .  (55) 


In  terms  of  these  parameters, 


dS'  =  Vg  do i  do%  ,  (56) 

with  g  the  first  fundamental  form  of  differential  geometry  for  S  in  terms  of  the 
parameters  o. 

In  Bleistein  [1987a],  the  method  of  stationary  phase  was  applied  to  the  integral 
(53)  in  the  variables  £  and  o.  The  conditions  of  stationarity  have  the  following 
interpretation.  First,  the  point  x'  on  S  is  such  that  Snell’s  law  is  satisfied  by  the  rays 
from  zT  and  xg  to  x\  A  particular  source/receiver  pair  is  picked  out  by  the  second 
condition  which  ties  the  triple,  x',  z,  and  zs  to  x.  This  condition  states  that  the  sum  of 
traveltime  gradients  from  z  to  xr  and  zt  and  the  same  sum  from  x  must  have  equal 
projections  on  the  observation  surface  defined  by  (46).  For  x  on  the  reflector,  S,  z  =z 
is  a  stationary  point,  with  Zj  and  z,  being  the  source  receiver  pair  whose  ray 
trajectories  satisfy  Snell’s  law  for  reflection  at  this  point.  The  condition,  /i=*=0, 
guarantees  that  there  is  only  one  such  point  for  sufficiently  smooth  surfaces,  S.  It  is 
assumed  that  the  surface  under  consideration  is  in  that  class.  Clearly,  then,  for  x  in  a 
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neighborhood  of  5,  there  is  still  only  one  stationary  point,  having  *'  =  x  as  its  limiting 
value  as  x  approaches  S. 

Let  us  suppose,  now,  that  x  is  near  S.  Introduce  a  neutralizer  function,  N(x\x), 
which  is  identically  equal  to  unity  in  some  disk,  say  |  x'  —  x  |  <  r^,  identically  equal  to 
zero  outside  some  larger  disk,  |  x'  —  x  j  >  r2  >  rlt  and  is  infinitely  differentiable 
everywhere.  Then  the  following  is  true. 

Lemma  2: 

Consider  the  integrals  obtained  from  (53)  by  introducing  N(x\x)  and 
N  (x',x)  =  1  —  N(x’,x)  as  multipliers  in  the  integrand.  Call  the  first  integral 
01t{x)  and  the  second  integral  /?jy(x).  Then  0s{x)  —  o(w~m),  where  m  is  limited 
only  by  the  smoothness  of  the  integrand. 

The  proof  is  given  in  Appendix  B. 

Now  consider  the  integral,  Ps(x).  This  integral  can  be  reduced  to  an  aperture- 
limited  Fourier  identity  operator  of  the  type  discussed  in  the  previous  section.  To  do 
so,  First  observe  [Beylkin,  1985]  that  for  x  near  x: 

u$(x,  x  ,  x„  xr)  =  u  V*  [r(x,  xr)  +  r(x,  xr)  •  (x'  -  x)  +  0  |  x  -  x  | 2]  .  (57) 

L  J  J  J 

That  is,  the  phase  is  nearly  of  the  form  lr(x'  —  x),  with 

k  =  wV,  |  r(x,  x,)  +  t(x,  xr)  .  (58) 

In  fact, 

Lemma  3: 

In  some  neighborhood  of  x  there  exists  a  change  of  variables  from  (u>,£)  to 
k  =  {kltk2,k3)  with  nonvanishing  Jacobian,  having  the  following  properties: 

k  =  w  V2  r(x,  xs)  +  r(x,  xr)  +0  |x'-x|  ;  (59) 

and 

=  uj2  H(x,()  =  [h(x,()  +  O  |x'-x|  .  (60) 

The  proof  is  given  in  Appendix  C. 
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It  is  interesting  to  note  that  the  regularity  of  the  transformation  does  not  depend 
on  cj,  but  only  on  the  spatial  variables.  Furthermore, 

Lemma  4: 

A 

k  —  k/\  k  |  is  a  function  of  £  sgnu  with  nonvanishing  Jacobian.  Conversely,  for 

each  choice  of  sgn  a/,  £  is  a  function  of  k. 

A 

That  is,  if  we  were  to  pick  two  parameters,  such  as  the  first  two  components  of  k  or  the 
polar  angles  of  k  with  respect  to  some  set  of  axes,  then  these  variables  are  functions  of 
£  (and  not  u)  with  nonvanishing  Jacobian  wherever  H(x,£)  is  nonvanishing,  except  for 
the  sign  of  u>.  This  follows  from  the  proof  of  Lemma  3. 

This  result  has  an  important  implication  with  regard  to  the  aperture  of  k  values  in 
the  domain  of  integration  after  transformation  to  these  variables.  The  directions  of  the 
k- vectors  in  the  domain  of  integration  is  solely  a  function  of  the  source/receiver 
configuration.  The  approximate  form  of  the  transformation  in  (59)  suggests  further 
that  this  angular  aperture  is  almost  completely  determined  by  the  sum  of  the  gradients 
of  the  travel  times  at  x.  Equivalently,  this  is  the  sum  of  the  tangents  to  the  rays  from 
the  source  and  receiver.  In  Beylkin,  Oristaglio  and  Miller  [1985]  and  Miller,  Oristaglio 
and  Beylkin  [1987],  extensive  use  is  made  of  this  fact.  For  our  own  purposes,  here,  the 
significance  of  this  observation  is  what  it  implies  about  the  integrand  of  (53)  after 
transformation.  Namely,  the  integral  takes  the  form 


/?(*)  =  JT  f  dakF(u/(k))  f  dS'  /(*',  k)  exp 
Sr1  J  Dk  ’'p 


ik-(z'  —  z) 


(61) 


In  this  equation, 


/(x',x,i)  =  -R(z’,zs) 


A(z',zs)A{z',zr)h-V’ 

t(z’,zs)  +  t{z\zt) 

1  0  1 

A{Z,Z,)A(Z,  Zr) 

Vr(x,  zt)  +  Vr(x,  xr) 

H{z,z',{) 

The  integral  (61)  is  exactly  of  the  form  of  the  integral  (44),  except  that  in  the 
application  here,  the  explicit  dependence  on  the  length  scale,  L ,  has  been  neglected. 
The  condition  that  the  wave  number  be  large,  in  order  to  apply  the  results  of  the 
previous  section,  first  requires  that  V[r(z,z,)  ■+■  r(x^xr)]  be  nonzero.  (Note  that  this 
gradient  vanishes  when  the  rays  from  z,  and  zr  are  tangent  at  x.  That  is,  they  are 
both  part  of  a  single  ray  from  z,  to  xr.  Such  ray  paths  would  naturally  arise  in 
transmission  tomography.  This  theory  does  not  apply  to  those  applications.  On  the 
other  hand,  diffraction  tomography  falls  in  the  range  of  applications  of  this  theory.) 
That  this  sum  of  gradients  does  not  vanish  is  assured  by  the  assumption  that  h(z,£)  be 
nonzero,  because  this  vector  is  the  first  row  of  the  determinant  in  (50).  Once  the 
minimum  magnitude  of  this  vector  is  determined,  and  the  natural  length  scale  L  for  the 
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integral  (53)  is  established,  the  burden  of  “large  wave  number”  is  put  on  the  frequency 
parameter.  In  seismic  applications,  frequencies  as  low  as  4  Hz  will  prove  to  be  high 
enough  for  the  corresponding  parameter,  A  =  KL,  to  be  large.  Thus,  the  theory  of  the 
previous  sections  applies  in  seismic  applications. 

It  follows  from  the  analysis  of  the  previous  section  that  the  output  of  (61)  is 
proportional  to  a  band  limited  singular  function  of  the  reflecting  surface.  The 
proportionality  factor  may  be  constructed  from  the  discussion  below  (44)  and  the 
definition  of  /  in  (52).  When  *  is  on  S,  this  factor  reduces  to  /  (x,x,Jfc).  However,  the 
distinguished  value  of  k  must  be  ±n  as  a  consequence  of  applying  stationarity  to  the 
derivatives  in  (14).  (This  condition  is  just  Snell’s  law.)  From  the  approximate  form  of 
k  in  (58) ,  it  is  apparent  that  k  must  point  downward ,  while  the  normal  to  5,  n,  in  the 
Kirchhoff  representation  of  U5  is  an  upward  normal.  Consequently,  in  (62) 


r(z',za)  +  t(x\Xt) 

=  — 

Vr(x,  z,)  +  Vr(x,  xT ) 

X  =x 

(63) 


Furthermore,  from  (60), 


0  =  •  (64) 

Thus,  when  the  distinguished  point,  x  on  S,  is  in  the  region  where  the  neutralizer 
functions  are  equal  to  unity, 


k)  =  i?(x,za)  .  (65) 

In  this  equation,  xs  is  fixed  by  the  stationarity  conditions  [(14)  and  the  change  of 
variables]  such  that  Snell’s  law  is  satisfied.  This  is  an  angularly  dependent  reflection 
coefficient.  It  is  a  function  of  cos©,  with  9  the  supplement  of  the  angle  between  n  and 
Vr(x,Xj).  At  the  stationary  point,  because  Snell’s  law  is  satisfied,  this  is  just  half  the 
angle  between  Vr(z,zs)  and  Vr(z,zr).  The  function,  cos©  is  determined  in  a  manner 
similar  to  the  way  cos#  and  coscf)  were  determined  in  the  previous  section,  using  the 
identity 


Vr(z,  za)  +  Vr(z,  zr)  =  2cos©  /  c(x) 


(66) 


which  cm  be  confirmed  by  squaring  both  sides.  See  Bleistein  [1987a]  for  details. 
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Application  to  Born- Approximate  Data 


In  Beylkin  [l985j,  the  Bom  approximation  was  used  to  represent  the  upward 
scattered  data.  The  inversion  of  such  data  for  determination  of  medium  parameters  is 
more  closely  akin  to  the  analysis  of  aperture-limited  Fourier  transforms  of  piecewise 
smooth  functions.  For  such  data,  we  already  know  from  the  previous  sections  what  the 
nature  of  the  aperture  limited  inversion  of  the  data  will  be.  For  Beylkin’s  inversion 
operator,  the  output  will  be  band  limited  step  functions  across  the  discontinuity 
surfaces  of  the  integrand  (usually,  a  perturbation  in  the  index  of  refraction,  or  some 
equivalent-  By  making  the  modifications  equivalent  to  changing  the  operator  /  to  the 
operator  I  of  the  previous  section,  the  output  will  again  be  singular  functions  of  the 
discontinuity  surfaces  scaled  by  the  jump  in  the  integrand. 

Thus,  whether  one  applies  the  operator  to  Bom  data  or  Kirchhoff  data,  the  result 
is  essentially  the  same,  a  weighted  singular  function.  The  only  issue  is  whether  one 
interprets  the  weight  in  terms  of  the  Bom  approximation  or  the  Kirchhoff 
approximation.  It  is  my  view  that  the  Kirchhoff  approximation  better  represents  the 
upward  scattered  field  in  most  applications.  One  reason  is  that  implicit  in  the  large 
wave  number  asymptotics  is  the  assumption  of  high  frequency  as  seen  above.  The 
Kirchhoff  approximation  is  a  high  frequency  approximation.  Furthermore,  in  many 
applications,  the  change  in  medium  parameters  across  reflectors,  or  the  angle  of 
incidence  and  reflection  with  respect  to  the  normal,  may  be  large  enough  to  call  into 
question  the  use  of  the  Bom  approximation. 

Summary 


It  has  been  shown  here  that  the  process  of  applying  a  proposed  inverse  scattering 
operator  to  model  data  is  asymptotically  equivalent  to  aperture-limited  Fourier 
inversion  with  a  dual  set  of  Fourier  variables  determined  from  the  travel  time  functions 
from  source  to  output  point  and  receiver  to  output  point.  The  inversion  operator 
proposed  here  is  a  modification  of  an  operator  proposed  by  Beylkin  [1985j.  Where  his 
inversion  would  lead  to  band  limited  step  functions,  this  one  leads  to  band  limited 
singular  functions.  When  the  bandwidth  is  of  sufficient  extent,  his  operator  leads  to  a 
Fourier  approximation  of  the  unknown  function.  When  the  bandwidth  is  not  sufficient 
for  this  application,  but  can  be  characterized  as  “large  wave  number”  in  some  sense, 
the  inversion  proposed  here  still  provides  a  map  of  the  discontinuity  surface(s) 
(reflectors)  of  the  unknown  function,  as  well  as  a  means  of  estimating  reflection 
strength. 

In  Bleistein  [1987a, b],  these  results  are  confirmed  by  applying  asymptotic  analysis 
directly  to  the  integral  (53).  Here,  it  was  shown  by  classical  asymptotic  methods  that, 
in  fact,  this  result  is  really  a  more  fundamental  phenomenon  of  Fourier  inversion  and 
the  results  for  inverse  scattering  follow  from  those  more  fundamental  results,  once  one 
identifies  the  Fourier  variables. 
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For  seismic  applications,  a  reflector  map  would  seem  to  be  a  desirable  output  of  an 
inversion  operator.  However,  Beylkin  and  his  associates  [in  the  above  cited  references] 
have  certainly  made  a  strong  case  for  a  more  complete  inversion  when  the  data  is 
available.  The  fact  that  the  operator  is  deduced  on  the  basis  of  a  high-frequency 
assumption  seems  not  to  degrade  the  quality  of  the  low  frequency  output  to  any 
significant  degree. 


CONCLUSIONS 


It  has  been  shown  here  by  classical  means  that  high  frequency  inversion  of 
scattered  data  is  equivalent  to  large  wave  number  aperture-limited  Fourier  inversion. 
The  type  of  output  one  obtains  from  the  inversion  operators  of  interest  can  be 
explained  in  terms  of  elementary  ideas  about  such  Fourier  inversion.  In  particular,  the 
transition  from  aperture-limited  step  function  output  to  aperture-limited  singular 
function  output  (across  the  steps)  follows  directly  from  the  asymptotic  analysis  of  large 
wave  number  aperture  limited  Fourier  transforms. 
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FIGURE  CAPTIONS 


Figure  1  The  singular  function  of  the  line,  2  =  0.  Also,  a  reference  “spike,”  to  be  used 
for  amplitude  comparisons. 

Figure  2  The  filter  function  applied  in  k-domain.  The  support  of  this  function  is  the 
domain  Dk. 

Figure  3  The  aperture  limited  inversion  back  to  the  spatial  domain  of  the  singular 
function  data  from  Figure  1,  restricted  to  the  domain  Dk  of  Figure  2. 
Reference  spike  in  the  comer  for  amplitude  comparison. 

Figure  4  An  alternative  filter  function  in  the  k-domain.  Dk  in  this  case  does  not 
include  the  direction  (0,1)  of  the  normal  to  the  support  of  the  singular 
function  in  Figure  1. 

Figure  5  Inversion  of  aperture  limited  data  for  the  singular  function  of  Figure  1  using 
the  domain  D k  of  Figure  4.  Reference  spike  demonstrates  that  the  scale  is 
the  same  as  in  Figures  3  and  5.  Inversion  is  effectively  zero  as  predicted  by 
the  theory. 

Figure  6  Same  output  as  Figure  5,  but  scaled  up  by  six  orders  of  magnitude. 

Figure  7  The  singular  function  of  a  circle. 

Figure  8  A  filter  containing  all  angles,  limited  magnitude  of  k.  The  domain  Dk  is  an 
annulus. 

Figure  9  Aperture  limited  version  of  the  singular  function  of  a  circle  using  the  domain 
Dk  of  Figure  8. 

Figure  10  An  alternative  aperture  to  apply  to  the  Fourier  transform  of  the  singular 
function  of  a  circle. 

Figure  11  Aperture  limited  inversion  of  Fourier  data  for  the  singular  function  of  a 
circle.  The  circle  is  properly  identified  only  for  the  region  in  which  the 
normal  to  the  circle  lies  in  the  angular  aperture  of  Dk. 
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APPENDIX  A 


Determination  of  the  determinant  and  signature  of  the  matrix  defined  by  (16)  will 
be  discussed  here.  The  elements  of  the  matrix  are  determined  from  (14)  and  (15). 
They  are  given  by 

d2$  -  d2y'  d2$  -  dy' 

doidoj  ^  doidoj  ’  80  do,  da,  ’ 


=  -  sin 0$>- ,  i,  j-  1,  2 

0<P  C/Oj  OO j 


32<5  _  32* 

-  =  —  p  •  (y  —  y ' )  -  =  —  cosff  <b  •  ( y  —  y ' ) 

a*2  *  y  j  ’  d0d<*  v  {V  y  } 


— —  = -p-(y-y')  ,  p  =  (cos<£,  sin0)  .  (Al) 

d<p 

With  no  loss  of  generality,  one  can  assume  for  this  calculation  that  ox  and  a2  are  each 
arclength  variables  in  the  principal  directions  at  the  stationary  point.  In  this  case,  the 
upper  left  two-by-two  matrix  of  elements  simplifies  to 


d2$ 
da,  da; 


=  -6ij  p-  Kj  ,  i,  y  =  1,2 


Here,  <5I;  is  the  Kronecker  delta,  equal  to  unity  for  » =  j,  zero,  otherwise.  Also,  in  this 
case,  the  first  derivatives  with  respect  to  a,-  in  (Al)  reduce  to  orthogonal  unit  vectors. 
In  addition,  the  vectors  0  and  $  are  orthogonal  unit  vectors,  co-planar  with  the  unit 
tangents  as  a  consequence  of  the  conditions  of  stationarity.  (Recall,  the  normal  to  Sy- 
is  colinear  or  anti-colinear  with  p.)  Therefore,  to  simplify  notation  below,  introduce  t /> 
as  follows: 


0--|SL-=cosV>,  ^•-^-=sin0  , 

do  i  do j 


dy' 

— * —  =  —  si 


=  —  sin0  ,  $ 


=  COS0  . 
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Direct  calculation  now  leads  to  the  result, 


det 


=  sin2# 


l-HxHiKX\y-y' 


l-/i2M3*2  |y  -  V'  i 


(A4) 


In  this  equation,  the  Hj’s  are  defined  by  equation  (20)  and  the  k}'s  are  the  principal 
curvatures  as  described  in  the  text.  This  result  can  be  seen  to  be  independent  of  rp. 
Thus,  for  the  purpose  of  determining  the  signature  of  $,y,  0  may  be  set  equal  to  zero. 
The  zeroes  at  0  =  0  or  n  are  due  to  the  singularities  of  the  coordinate  system  at  these 
points  and  can  be  disregarded  here. 

It  cam  be  seen  from  (A4)  that  the  sgn  $,y  can  only  change  when  y  moves  through 
a  center  of  principal  curvature  of  Sy-.  Thus,  for  the  purpose  of  determining  sgn  4>,y  for 
y  near  Sy-,  we  cam  take  y  on  Sy-\  that  is,  we  may  set  |  y  —  y  |  =0.  In  this  case,  we  find 
that 


det 


$«y  My 


=  sin 


A2  ■+•  —  1 


A2  +\fi2K2  —  sin2# 


(A5) 


Each  factor  here  has  two  roots  of  opposite  sign,  leading  to  the  conclusion  that  sgn 
$tJ  =  0.  This  completes  the  proof. 


APPENDIX  B 

Lemma  2  will  be  proven  in  this  appendix.  The  fourfold  integral  /?yy(i)  in  £  and  a 
is  of  the  form 

Pn{z)  =  Jd  G{*l)  exp  ^(ii)  j  dr^d^  dr]zdr)A  .  (Bl) 

In  this  equation,  the  four  variables,  17,  are  just  the  variables,  (£1,  renamed  to 

make  the  discussion  below  easier.  The  phase  is  just  the  phase  $  in  the  newly- 
named  variables,  while  G  identifies  the  full,  neutralized  amplitude  of  (54)  multiplied  by 
N  (x‘,x).  The  domain,  Dn,  is  just  the  pair  of  domains  in  £  and  a:  S^xSa.  The 
dependence  of  the  integrand  on  x  has  been  surpressed  because  it  is  not  important  in 
this  discussion. 

The  integrand  has  no  stationary  points  in  Dv  and  vanishes  C°°  smoothly  on  its 
boundary.  Therefore,  the  integrand  is  expanded  as  was  done  in  the  previous  section 
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eq.  (9)],  in  preparation  for  integration  by  parts: 


r  11  v  r  i  ci 

G(i7) exp  iu$i{ri)  -  —  V„*  - 2*  G(»j)  exp  i\^  +  Gt(  17)  exp  iA^ 

L  tu>  V$x  L  J  L  J 


Gi(n)  =  -  v„ 


"W  G(^) 


Substitute  this  identity  into  the  previous  equation  and  use  the  divergence  theorem  to 
replace  the  first  integral  over  Dn  by  an  integral  over  the  boundary  of  this  domain. 
Since  the  integrand  vanishes  on  the  boundary,  this  first  integral  is  zero  and 

Pn(x)  =  ~~  fD  Gi(i?)exp  [twfcifo)]  drjx  dri2  drj3driA  .  (B3) 

This  process  can  be  repeated  recursively,  as  long  as  the  integrand  has  derivatives  in 
Dv.  Each  time,  another  power  of  l/t'w  is  introduced  as  a  multiplier  of  the  integrand 
while  the  integral,  itself,  remains  finite.  Thus,  the  claim  of  the  lemma  is  confirmed. 


APPENDIX  C 


Lemma  3  will  be  proven  here.  For  brevity,  set 


'Ffc  *»,  *r)  =  *«)  +  t(x,  Xr)  ,  (Cl) 

and  write  the  Taylor  series  for  as 

*(*,  x„  x,)  =  V*(x,  x.,  x.)  +  £j  £  i  W  •?’*}.  (,•  -  ,)»  .  (C2) 

n  =  2  Ill'll  =  n 

In  this  equation, 


v  =  fa.**.**)  >  lli/l!  =  vx  +u2  +v3  , 
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v\  =  I/! !  I/j !  1/3 !  ,  (*'  -x)V  =  {x{  -XiY1  (12  -i2)‘/l  (*3 

ai"  =  dx\'  dx?  dxf  ,  (C3) 

with  the  Uj’s  being  nonnegative  integers. 

Now  define  the  vector  k  by 

kj  =  uipj  , 

n  =  d*  v  v1 

P;  a*y  +  ii„ii  I/! 

•'  n  =  2  lli/ll  =  n 

(C4) 

One  can  check  that  k-(x'  —  x)  equity  by  multiplying  p;  by  w(ajp  xy  and  summing  on  j  to 
obtain  the  series  in  (C2)  multiplied  by  u>.  Note  that  the  radius  of  convergence  of  each 
of  the  series  in  (C4)  is  the  same  as  for  the  series  in  (C2).  Furthermore,  (59)  follows  by 
direct  computation.  Thus,  in  some  neighborhood  of  x,  the  transformation  from  (w,0  ^ 
x  is  one-to-one.  Indeed,  with  hindsight,  choose  N(x',x)  to  be  small  enough  that  H(x,£) 
is  nonvanishing  on  the  support  of  N(x’,x). 

This  completes  the  proof. 
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